Extraction of Spin-Dependent Parton Densities and Their Uncertainties 
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We discuss techniques and results for the extraction of the nucleon's spin-dependent parton dis- 
tributions and their uncertainties from data for polarized deep-inelastic lepton-nucleon and proton- 
proton scattering by means of a global QCD analysis. Computational methods are described that 
significantly increase the speed of the required calculations to a level that allows to perform the full 
analysis consistently at next-to-leading order accuracy. We examine how the various data sets help 
to constrain difi^erent aspects of the quark, anti-quark, and gluon helicity distributions. Uncertainty 
estimates are performed using both the Lagrange multiplier and the Hessian approaches. We use 
the extracted parton distribution functions and their estimated uncertainties to predict spin asym- 
metries for high-transverse momentum pion and jet production in polarized proton-proton collisions 
at 500 GeV center-of-mass system energy at BNL-RHIC, as well as for W boson production. 
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I. INTRODUCTION 



The last twenty years have witnessed remarkable im- 
provements in the sophistication and precision of meth- 
ods devoted to the extraction of parton density and frag- 
mentation functions from experimental hard scattering 
data. These distributions are essential ingredients for 
any phenomenological study of hard scattering processes 
involving identified hadrons in the initial and final-state, 
respectively. Their precise knowledge is not only critical 
for testing the very successful framework of perturbative 
Quantum Chromodynamics (pQCD) but, in more gen- 
eral terms, also for quantifying uncertainties for preci- 
sion studies of the Standard Model and searches of "new 
physics" at high energy accelerators like the CERN-LHC. 
At the same time, parton densities and fragmentation 
functions give fundamental insights into nucleon struc- 
ture and the hadronization mechanism. 

With the gain in precision and refinements of analyses, 
modern parton distribution functions (PDFs) have often 
revealed intriguing aspects of hadronic structure, such 
as the sizable breaking of isospin symmetry in the light 
sea quark sector, suggestions of differences between the 
strange quark and anti-quark distributions, the steep rise 
of the distributions at small momentum fractions, and an 
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interesting pattern of modifications of the distributions 
in nuclei, to name just a few. Certainly one of the most 
striking results is the unexpectedly small fraction, about 
a quarter, of the proton's spin that can be attributed to 
the intrinsic angular momenta of quarks and anti-quarks. 
This finding, famously dubbed "proton spin crisis" , has 
triggered a flurry of experimental and theoretical activ- 
ity aiming at clarifying the contributions of gluons and 
orbital angular momenta of partons to the spin of the 
proton 1]. 

The only way to effectively deconvolute the experimen- 
tal information on PDFs, which in its raw form is smeared 
over the light-cone momentum fraction x, summed over 
many different partonic subprocesses, and taken at dif- 
ferent hard scales Q for each data point, is a "global 
QCD analysis" . It treats all available probes simulta- 
neously, in order to extract the set of universal PDFs 
that yields the optimal theoretical description of the com- 
bined data. For the case of polarized PDFs, the avail- 
able world-data are from polarized deep- inelastic scat- 

SM, [l^ [H [H, [il , semi- 
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tering (DIS) [3, ^, 1, d 

inclusive DIS (SIDIS) [io|, Ijj, ll5|, [lelLphoto- and elec- 
troproduction of hadrons and charm l3, [H, [H, [20l. [2ll 
and proton-proton {pp) collisions at BNL-RHIC [22, [2j, 
HSj H^, [13 • The different data sets are complemen- 
tary in the sense that they probe different aspects of the 
helicity dependent PDFs. Fully inclusive DIS data from 
the many different experiments are pivotal in precisely 
determining the sums of quark and anti-quark distribu- 
tions, SIDIS data help to tell different quark flavors and 
quark and anti-quarks apart, and RHIC pp data give a 
first direct constraint on gluon polarization. 
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A global QCD analysis of nucleon spin structure at full 
next-to-leading (NLO) accuracy was completed recently 
[28l |. The present paper gives in large part a more de- 
tailed account of the methods and results of [28^1 . It also 
addresses the issue of the uncertainties of the PDFs in a 
more detailed and comprehensive wa y. A s customary for 
recent unpolarized PDF analyses [29l.l30|| , we provide sets 
of polarized PDFs associated with displacements in the 
PDF parameter space in the vicinity of the best fit which 
greatly facilitate the propagation of PDF uncertainties to 
any observable of interest. 

As a new feature over all previous fits based only on 
DIS m, m, UlLor combined DIS and SIDIS 

data [3y|, the analysis 2^ included for the first time 
also results from polarized pp scattering at RHIC in 
a NLO framework. It benefitted significantly from an 
improved knowledge of parton-to-hadron fragmentation 
functions which are an essential non-perturbative in- 
put for the theoretical description of all processes with 
identified hadrons in the final state, such as SIDIS. For 
the first time, these fragmentation functions provide a 
good description of identified hadron yields in the entire 
kinematic regime relevant for the analysis of polarized 
SIDIS and pp data [S^l • For the time being, hadron and 
charm pro duction data from fixed-target experiments 
[TtI . [TsI. Il9l [20I [2]| . which constrain the gluon polariza- 
tion at momentum fractions around a: ~ 0.1, were not in- 
cluded in the analysis [1^ since NLO calculations of the 
relevant cross sections are not yet complete [l^ . We will 
provide a comparison to these data in order to demon- 
strate their consistency with the results of the global fit. 
RHIC data for charged pion spin asymmetries [131 are 
also not taken into account as they are still preliminary 
and statistically not as significant as the neutral pion or 
jet data [H, Hl, [H, [H, [2i|. With sufficient statistics, 
however, they can provide an important constraint on 
the gluon polarization as will be shown below. 

The use of parton distributions in predictions for Stan- 
dard Model benchmark processes, e.g., as "luminosity 
candles" at the LHC, or in understanding fundamental 
properties of a nucleon like its spin, not only requires a 
careful extraction of PDFs from data but also a proper 
assessment of their uncertainties and how they propagate 
to other observables of interest. In spite of a great deal 
of activity and many significant achievements for both 
parton distribution and fragmentation functions (FFs), 
this has shown to be a rather formidable task in prac- 
tice [H, HO, [13, [M Hi- The specific challenge of a 
global QCD analysis is to incorporate a large body of 
data from many experiments with diverse characteris- 
tics and errors. The complications are compounded by 
uncertainties inherent to the theoretical framework used 
to describe the data, which are notoriously difficult to 
quantify. Examples are the choice of the factorization 
scale, the functional form used to parameterize the PDFs, 
or unavoidable approximations and assumptions limiting 
the parameter space. 

Several complementary strategies have been devised 



and implemented to estimate uncertainties of PDFs and 
FFs [U, [H, [4^ . In general, one starts with introducing 
an effective function that combines all phenomenologi- 
cal inputs to the analysis as a quantitative measure of the 
goodness of the global fit. Minimizing this function 
yields the optimal set or "best fit" of parameters in the 
multi-dimensional space defining the PDFs. The most 
common method to determine the range of uncertainties 
is to study the dependence of near its global minimum 
based on a Taylor expansion and keeping only the leading 
term as characterized by the error matrix or its inverse, 
the Hessian matrix . This assumes a quadratic form 
in the displacements of all parameters from their opti- 
mum values. The Hessian also determines the uncertain- 
ties of any other physical observable O, provided that the 
dependence of O on the fit parameters is approximately 
linear around the minimum. Both assumptions are not 
necessarily adequate in the complex global analysis en- 
vironment, and their range of applicability needs to be 
carefully scrutinized. 

The more robust method of Lagrange multipliers [i^l 
circumvents all these shortcomings and is free of assump- 
tions concerning the functional dependence of on the 
fit parameters. The idea is to explore directly how the fit 
to data deteriorates if one enforces certain values for an 
observable O away from its best fit value. In practice, one 
performs a series of constrained fits in which x^ is mini- 
mized for particular values of O, in order to map out the 
parametric relationship between and O. The method 
is straightforward to implement and can be applied to 
any combination of physical observables or even to fit pa- 
rameters themselves. We will pursue and compare both 
methods, Hessian and Lagrange multiplier, to estimate 
uncertainties for the shape and truncated first moments 
of helicity-dependent PDFs using the analysis presented 
in [2^ as the starting point. The Lagrange multiplier 
technique will provide the necessary benchmarks for test- 
ing the accuracy of approximations within the Hessian 
method. We note that alternative approaches recently 
proposed in the literature include studies of uncertain- 
ties based on neural networks or large samples of PDFs 
generated with Monte-Carlo methods [ioj . 

In any case, all of these methods require an extensive 
number of calculations and minimizations of the effective 
X^ function, in order to explore the very complex and en- 
tangled sensitivity of the data to variations of the param- 
eters describing the PDFs. This calls for new and more 
efficient calculational tools to include all observables used 
in the global fit consistently at NLO accuracy without re- 
sorting to potentially unreliable approximations. In par- 
ticular, numerical computations of NLO cross sections in 
hadron-hadron scattering are prone to being very time 
consuming. 

In order to deal with this problem, Ref. [2^ employed 
a method based on the Mellin transform technique pro- 
posed in Refs. [11, [1^, which allows to speed up the rele- 
vant NLO computations to a level that they can be incor- 
porated in the global analysis. An important improve- 
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ment of this Mellin transform method was the implemen- 
tation of a novel Monte-Carlo sampling technique. This 
new computational strategy proved already very useful 
for the analysis of single-inclusive observables in pp scat- 
tering that are presently relevant at RHIC. However, it is 
completely general and becomes especially powerful when 
less inclusive observables like two-particle correlations in 
hadronic collisions need to be incorporated in the global 
QCD fits, as will soon be the case. In the present pa- 
per we will describe the Mellin transform technique and 
its improvement in detail. We note, that the fast and 
efficient Mellin technique for incorporating NLO pp pro- 
cesses is, of course, not restricted to analyzing helicity- 
dependent PDFs, but could equally find important ap- 
plications for QCD processes at the LHC. 

The paper is organized as follows: in the next Sec- 
tion we describe all technical details of a global PDF 
analysis. We first discuss the function and the un- 
derlying ideas of the Hessian and Lagrange multiplier 
approaches for estimating PDF uncertainties. We next 
describe in detail our Mellin moment and Monte-Carlo 
sampling techniques as implemented for fast evaluations 
of NLO pp cross sections in our global QCD analysis. 
In Sec. HI we apply all techniques to the global anal- 
ysis of helicity-dependent PDFs [1^. We discuss the 
results for the best fit and its uncertainties. We argue 
that the range of applicability of the Hessian method is 
limited to estimating uncertainties of helicity-dependent 
PDFs consistent with only small departures from the best 
global fit, corresponding to Ax^ ~ 1. We present sets of 
polarized PDFs associated with displacements along the 
eigenvector directions of the Hessian matrix and result- 
ing in Ax^ = 1, which characterize the PDF parameter 
space in the vicinity of the global minimum in a process- 
independent way. We also explore the impact of the indi- 
vidual data sets on the results and uncertainties obtained 
for helicity-dependent PDFs in the global fit. In Sec. IV 
we study the potential of upcoming measurements at 
RHIC a.t y/S ^ 500 GeV center-of-mass system (c.m.s.) 
energy for further constraining the polarized PDFs. We 
focus on predictions for single-inclusive pion and jet pro- 
duction, and on W boson single-spin asymmetries. We 
conclude in Sec. V. 



A. The effective function 

Global QCD extractions of PDFs [H, S M, S M, 
l47j or FFs [s^ are implemented around an effective x^ 
function that quantifies the goodness of the fit to data for 
a given set of theoretical parameters {fli}, i — I, . . . , A^par 
that determines the PDFs or FFs at some initial scale /io- 
The simplest function, convenient for the search for 
optimum PDFs by minimization, is usually taken as 



ri=l j = l 



D,-T,{{a,}) 



SDi 



(1) 



where iVoxp counts the individual experimental data sets 

(n) 

and -^^ata corresponding number of data points in 
each set. Each data value Dj is compared to the corre- 
sponding theoretical estimate Tj, which depends in gen- 
eral non-linearly on the TVpar parameters {a-i}, weighted 
with the estimated uncertainties combined in SDj. In 
Eq. ^ LOj is a special weighting factor for each data 
point with default value one. It can be set to zero if 
a certain data point is to be removed from the analysis 
due to some physics considerations. For instance, such 
cuts are routinely introduced in a global fit to remove 
kinematical regions where the framework of perturbative 
QCD used to compute Tj{{ai\) is known to be not ade- 
quate for describing the available data. The simple form 
([T]) for x^ is appropriate only in the ideal case of data 
sets with uncorrelated errors, and (5Z?J is then given by 
statistical and point-to-point systematic errors added in 
quadrature. 

For most experiments, additional information on the 
fully correlated normalization uncertainty 5Mn can be 
found, i.e., on a systematic shift common to the entire 
data set. Equation (jT]) is straightforwardly extended to 
account for such normalization uncertainties: 



J\fnDj-Tj{{a,}) 
6D, 



(2) 



II. TECHNIQUES FOR NLO GLOBAL PDF 
ANALYSES 



In this Section, we will describe all techniques we use 
for the global analysis of polarized PDFs. The first 
two Subsections discuss the x^ function and the various 
methods for the analysis of PDF uncertainties. Much 
of the discussion here will follow the pioneering work in 
Refs. m El, El- We then lay out the details of our 
Mellin moment and Monte-Carlo sampling techniques. 



Here, Afn are the normalization factors, which can be 
either fitted along with the {ai}, or even determined 
analytically in each step of the minimization by solving 

There are several equivalent methods of extending fur- 
ther the simple function in Eq. ^ in the presence of 
sources of correlated systematic errors for a data 
point Dj in data set n ^42, i46] ■ The numerically most 
efficient method treats the correlated systematic errors 
analytically in the optimization procedure like for the 
global normalization uncertainties discussed in ^ . This 
avoids the construction and inversion of large covariance 
matrices used in the conventional approach. The result- 
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ing function has the form (assuming LUj = 1 for sim- 




(3) 



where 



N 



Bk 



E 



SD 



(n) 



AT 



(4) 
(5) 



Here, (JD^"^)^ is the quadratic sum of the statistical and 
uncorrelated systematic errors, and f3kj specifies the fcth 
correlated systematic error of data point Dj . 

We note that most global analyses of unpolarized PDFs 
now start to include correlated systematic errors when- 
ever this information is available from experiment. This 
is of much importance, as very precise PDF uncertainty 
studies are a key ingredient for reliably estimating new 
physics signals and standard model backgrounds at the 
Tevatron or the LHC. For the time being, and keeping in 
mind that our analysis is the first fully global one of this 
kind for polarized PDFs, we stick to an effective func- 
tion based on the simplest expression in Eq. ([1]). Also, 
most of the relevant experiments do not publish the full 
information on correlated systematic errors. Whenever a 
global normalization uncertainty is provided, we have ex- 
plored the possibility of normalization shifts to improve 
the global fit by minimizing according to Eq. ^ . We 
have not found any significant improvements of the fit 
from this. Since data sets are continuously evolving and 
more and more precise information becomes available, 
the proper inclusion of correlated systematic errors is cer- 
tainly needed in future global analyses of helicity PDFs. 

As mentioned in the Introduction, there are consider- 
able complications when statistical methods are applied 
to global QCD analyses based on a large body of di- 
verse data and a theoretical model with many parameters 
{fli}. In particular, the statistical value of the definitions 
given in Eqs. (IT|)-(I^ has been under considerable de- 
bate [2i,|33,liiljS,lSlii,|43, since both the theoretical 
[Tjdai})] and the experimental inputs [Dj, 6Dj, . . .] are 
far from being ideally suited for statistical analysis. For 
instance, uncertainties inherent to the theoretical frame- 
work used to describe the data are notoriously difficult to 
quantify and usually correlated. In addition, it is often 
the case that even in a "best fit" to data, the values of 
per data point for individual experiments vary con- 
siderably around unity, sometimes by much more than 

the expected amount v'^/n'^^^^- 



drawn on a tolerable range of uncertainty from a certain 
increase A^^, must be considered keeping in mind this 
complex situation. 



B. Uncertainty estimates: Hessian and Lagrange 
multiplier methods 

An important objective is to estimate uncertainties of 
the spin-dependent PDFs obtained from the optimiza- 
tion. To this end one must study the behavior of the effec- 
tive x^({ai}) function in the neighborhood of the global 
minimum Xq = X^({*^i })j using reliable statistical meth- 
ods rather than some subjective tuning of selected pa- 
rameters. Basically two complementary tools have been 
put forward, refined to deal with the complexity of global 
PDF analyses, and applied to characterize uncertainties 
in the extraction of unpolarized PDFs. We pursue and 
compare both, the Hessian [i^ and the Lagrange multi- 
plier [i^ methods, in our uncertainty estimates. We only 
briefly recall the main elements of the two approaches 
and highlight potential problems and limitations. For a 
detailed account we refer the reader to Refs. [HI, [H, El] ■ 

In the more standard Hessian approach, the explo- 
ration of the uncertainties associated with the fit is imple- 
mented through a Taylor expansion of X^{{ai}) around 
the global minimum Xo({'*i})- Keeping only the lead- 
ing quadratic terms, the increase Ax^ can be written in 
terms of the Hessian matrix 



1 d\ 



2, ,2 



2 dy^dyj 



Ax' = x'({a.}) - xlM}) = E H^^y^y^ 



(6) 



(7) 



where yi = ai — and the derivatives in Eq. ([6]) are 
taken at the minimum. 

Global QCD fits are usually characterized by very dis- 
parate uncertainties in different directions of the multi- 
parameter space, so that standard methods to evaluate 
Hij by finite difference, as implemented in, e.g., the Mi- 
NUIT package tend to be numerically unstable and 
hence unreliable. To overcome such difficulties, a new 
iterative algorithm to compute the derivatives in ^ was 
devised in Ref. [ilj and subsequently used in global anal- 
yses of unpolarized PDFs. We apply this improved Hes- 
sian method also in our studies. 

Instead of working in the parameter basis {oi}, the 
Hessian Hij is expressed in terms of its A^par eigenvectors 

v'l'^^ and eigenvalues Sk- The displacements yi in Eqs. ([6]) 
and ([7]) are replaced by a new set of parameters {2^} 
defined by [4l|, ]^ 



Therefore, conclusions 



(8) 
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The {zi} are appropriately normalized by scale factors 
Sj oc \/l/ej such that surfaces of constant turn into 
hyper-spheres in {zi} space, with the distance from the 
minimum given by 



(9) 



Large eigenvalues Sk correspond to steep directions in 
which changes rapidly and the parameters are tightly 
constrained by the data, while small eigenvalues belong 
to directions where the parameters are only weakly de- 
termined. 

Within the eigenvector representation, it is convenient 
to construct 2iVpar eigenvector basis sets of PDFs which 
greatly facilitate the propagation of PDF uncertainties 
to arbitrary observables O [43|]. These basis sets are 
defined in {zi} space by 



(10) 



and hence correspond to positive and negative displace- 
ments along each of the eigenvector directions by the 
amount T = ^ Ax^ still tolerated for an acceptable 
global fit. To estimate the error AO on a quantity O 
away from its best fit estimate 0(5'°) it is only necessary 
to evaluate O for each of the 27Vpar sets S'^ [43|, i.e., 



AO = 



k=l 



1/2 



(11) 



One must keep in mind that the propagation of PDF 
uncertainties in the Hessian method has been derived 
under the assumption that a first order, linear approxi- 
mation is adequate. Of course, due to the complicated 
nature of a global fit, deviations, also from the simple 
quadratic behavior in Eq. ([7]), are inevitable, and error 
estimates based on the Hessian method are not necessar- 
ily always accurate. 

A strategy based on Lagrange multipliers avoids all 
these assumptions and probes the uncertainties on any 
quantity 0({ai}) much more directly. The result is a 
parametric relationship between one or more observables 
Oj({ai}) and the effective function used to determine 
the goodness of the global fit. Its application to global 
QCD analyses was proposed in Ref. [i^l not only to es- 
timate uncertainties of observables depending on PDFs 
but also to test the range of applicability of the Hessian 
approach described above. 

In practice, the method is implemented by minimizing 
a function 

*({«J,{A,}) = x^{{a,})+Y,^,0,{{a.,}) (12) 

3 

with respect to the set of PDF parameters {ci} for fixed 
values of the Lagrange multipliers {A^}. Each multiplier 
is related to one specific observable Oj, and the choice 



\j = corresponds to the best fit 5'". By repeating this 
minimization procedure for many values of \j one can 
map out precisely how the fit to data deteriorates as the 
expectation for the observable Oj is forced to change. 
Here, contrary to the Hessian method, no assumptions 
are made regarding the dependence of or the observ- 
able Oj on the parameters {a^} of the fit. The Lagrange 
multiplier method also generates a large set of sample 
PDFs along the curve of maximum variation of the ob- 
servable(s) used in Eq. p^ . 

In principle, the Lagrange multiplier method is supe- 
rior to the Hessian approach for reliably estimating un- 
certainties. For a given Ax^ it finds the largest range of 
Oj{{ai}) allowed by the data used in the global fit and 
the theoretical ansatz, independent of the approxima- 
tions involved in the Hessian method. The entire param- 
eter space {ai\ is explored in minimizing Eq. p2p . not 
necessarily limited to the vicinity of the best fit {a^}. In 
practice, a large number of global fits is required to map 
out all x^ profiles of interest. Thanks to our novel Monte 
Carlo sampling technique to be described below, this is 
no longer a serious limitation computationally. The only 
drawback to this method is that it requires access to the 
full machinery of global fitting to estimate uncertainties 
of a given observable of interest, contrary to the Hessian 
method for which the eigenvector PDF sets make it 
very simple to propagate PDF uncertainties to arbitrary 
observables. 

In this paper, Lagrange multipliers will mainly be the 
method of our choice for estimating uncertainties and to 
monitor to what extent the approximations involved in 
the Hessian approach are justified. 



C. Computational technique: The Mellin method 

As is evident from the previous Subsections, the ex- 
traction of PDFs in a global QCD analysis of a large body 
of data at NLO accuracy requires an extensive number of 
time consuming computations of the corresponding ob- 
servables in each step of the x^ minimization procedure. 
The large number of parameters specifying the functional 
form of the PDFs in the fit, typically of 0(20), and the 
need for a proper assessment of their uncertainties, add 
to this. Assuming a representative analysis with about 
500 data points, 5000 iterations to reach the optimum 
set of parameters, and a modest computing time of 1 
second per cross section calculation at NLO, one easily 
realizes that computational improvements are very much 
in demand. 

We stress at this point that approximating NLO cor- 
rections by parameterizing them by a -fC-factor, K = 
^NLO^i^LO^ which is a possible strategy in the spin- 
averaged case in order to speed up the analysis [1^ H^, 
m, So, 111,113) is not a viable option in the polarized case. 
In the latter, the parton distributions as well as the hard 
scattering cross sections may have nodes and change sign 
within the kinematic region of interest. As a result, differ- 
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ent partonic subprocess contributions can have very dif- 
ferent NLO corrections that are never well approximated 
by a single X-factor. This problem is even more promi- 
nent at the present time when the spin-dependent par- 
ton distributions, in particular the giuon distribution, arc 
still poorly known. Also, from a more theoretical point 
of view, the NLO corrections are expected to decrease 
the factorization/renormalization scale dependence of the 
calculation. In the ideal case that the NLO cross section 
has relatively little scale dependence, this would imply 
that the i^T-factor inherits the full LO scale dependence 
and thus cannot serve well as a proxy for the NLO cor- 
rections. We thus do not really see a useful workaround 
that would avoid inclusion of the full NLO calculation in 
the global analysis. 

The required numerical calculations also involve the 
scale evolution of the PDFs from some initial scale ~ 
0(1 GeV) to each of the scales relevant for the data 
points. The evolution is governed by a well-known set of 
integro-differential equations [H, [SO] that can be solved 
analytically after an integral transform from Bjorken x 
space to complex Mellin iV-moment space. The Mellin 
transform of a generic function tp depending on the light- 
cone momentum fraction x is defined as 



$(iV) = / x'^-'^(p{x)dx, 
Jo 



and its inverse reads 



ip{x) 



1 
27ri 



a;-^$(iV) dN. 



(13) 



(14) 



Here Cn denotes a suitable contour in the complex N 
plane that has an imaginary part ranging from — cx3 to 
-l-oo and that intersects the real axis to the right of the 
rightmost pole of <i>(iV). In practice, it is beneficial to 
choose the contour to be bent at an angle < 7r/2 to- 
wards the negative real-A^ axis [5l|. The integration in 
()14p can then be very efficiently performed numerically 
by choosing the values of N as the supports for a Gaus- 
sian integration. 

The transformation has the welcome property 

that convolutions factorize into ordinary products, which 
greatly simplifies calculations based on Mellin moments. 
It can be carried out analytically not only for the splitting 
functions governing the scale evolution of the PDFs but 
even for the partonic hard scattering cross sections for 
both inclusive and semi-inclusive DIS. The latter case re- 
quires straightforward extensions of Eqs. and C5|to 
double transformations as was discussed in Ref . [44, . 
This reflects the fact that SIDIS depends on two scaling 
variables, the momentum fractions x and z taken by the 
struck parton from the parent nucleon and by the final- 
state hadron from the scattered parton, respectively. The 
usefulness of the Mellin technique was demonstrated in 
practice in a global analysis of helicity-dependent PDFs 
using all available DIS and SIDIS data |36| . 

However, the inclusion of observables in hadron-hadron 
collisions or in less inclusive reactions in lepton-hadron 



scattering, which are crucial for determining, e.g., the 
gluon distribution, require a more elaborate framework. 
They involve multiple convolution integrals between one 
or more PDFs, partonic cross sections, and, depending on 
the process, fragmentation functions. More importantly, 
they typically depend on several kinematic variables, so 
that there is no direct way of taking Mellin moments of 
the cross section under which it would become a simple 
product of Mellin moments of PDFs and partonic cross 
sections. An example is single-inclusive pion production 
in proton-proton collisions, pp —^ nX, at measured trans- 
verse momentum px and rapidity rj of the pion. While 
XT — '2pT I ^/S is the typical scaling variable of the pro- 
cess, taking Mellin moments of the cross section in x\ 
does not lead to a simple expression involving the mo- 
ments of the PDFs. An efficient computational scheme 
that allows to circumvent this complication has been de- 
vised in Ref. f43,'45'|. Its main feature is to use the inverse 
Mellin transforms of the PDFs in the theoretical expres- 
sion, which makes it possible to store all numerically time 
consuming calculations involving the lengthy and compli- 
cated expressions for the underlying hard scattering cross 
sections on large "look-up tables" or "grids" in Mellin 
moment space. Here, we review the technique of [4^ 1 and 
describe a method to compute these grids very efficiently 
using Monte-Carlo sampling techniques. The latter al- 
lows to generalize the Mellin moment technique beyond 
the single inclusive processes considered in Ref. ^1] . 

A typical observable of interest, the spin-dependent 
cross section for pp — > -kX at RHIC, has the following 
general factorized structure: 



AaL = -[AaL. (++)-Aa|,. (+-)] 
^Y.\ A/,;(xi)A/,(x2)i?fc(z) 



(15) 



X dlS.Uij^kx {xii X2, z)S dxidx2dz . (16) 



Here we have suppressed the dependence on kinematic 
variables other than momentum fractions, and on the fac- 
torization and renormalization scales. In (jl6p . Afi and 
Dk denote the spin-dependent parton distribution and 
spin-averaged fragmentation functions for flavor i and k, 
respectively, and dAaij^kx the relevant partonic hard 
scattering cross sections. S represents a "measurement 
function" that defines details of the observable Act |bin 
in each bin, such as the kinematical ranges explored 
and the relevant experimental cuts. The symbols ± in 
Eq. (fTS]) denote the helicity combinations of the colliding 
longitudinally polarized protons defining the cross sec- 
tion Act Ibin- 

Following the strategy developed in Ref. we re- 
place the PDFs in Eq. by their representations as 
Mellin inverses, Eq. and subsequently interchange 
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integrations to obtain 

ijk •^'^M 

xdAa^jk{N,M)dNdM , (17) 

where 

i:,,(Ar,Af)EE A/,(Ar)A/,(M) (18) 

and 

dAa^jkiN,M) EE J dxidx2dzx^^x^^'Dkiz) 

X dAaij^kx{xi,X2,z)S . (19) 

Here and in the following, we assume that only the spin- 
dependent PDFs are subject to the global analysis and 
that the fragmentation functions are known. 

In the next step, one can now calculate the quantities 
dA<Tijk{N, M) in Eq. (fT9|l prior to the actual fit, as they 
do not depend on the PDFs, and store them in large look- 
up tables in the Mellin variables M and N along the con- 
tours Cm,n- In practice, it is convenient to choose M and 
N on the supporting points of a Gaussian integration, see 
Ref. [31 for details. As can be seen from (|19p. this effec- 
tively amounts to computing the NLO cross sections for 
all partonic subprocesses using complex "dummy PDFs" 
Xi^x^^^ ■ Even after making use of all symmetry rela- 
tions, e.g., among different subprocesses ij — > kX, set- 
ting up all look-up tables for a typical N x M grid size 
of 64 X 64 is a rather time-consuming step in practice, al- 
though it has to be done only once. This is where our new 
Monte-Carlo sampling technique has considerable advan- 
tages over a "brute-force" computation. 

First, we recall that a Montc-Carlo algorithm reduces 
the multiple integrations in (jisp to a finite sum over ran- 
dom "events" n, 

(20) 

n=l 

with weights w*^"^ and scaled by the number of itera- 
tions, K, used to generate large samples of / events. The 
weight includes the dependence on the parton distribu- 
tion and fragmentation functions and the hard scattering 
cross sections for each event. In Eq. (f20l) we keep the 
measurement function S separate from as it is usu- 
ally implemented only after an event has been generated. 

The most efficient strategy to compute all the re- 
quired dAaijk{N, M) in Eq. (fT9|) is to choose a suit- 
able set of trial PDFs Afi and to perform a single high- 
statistics Monte-Carlo integration calculation of the cross 
section in Eq. (fT6|) . During the calculation one stores, 

(n) 

for each event n, the momentum fractions x\ 2 and the 

corresponding weights w^^l for all of the subprocesses 
ij — > kX. Renormalizing each weight by 

-S'^-S^i? (21) 



with 

L[f^AMx^r^)Af,{x^-^), (22) 

removes its dependence on the assumed set of PDFs. The 
dAdijk {N , M) are then straightforwardly obtained as 

n=l 

(23) 

In other words, knowledge of the set a;^"-* , Xj"'', w^"^'' 
which corresponds to a profile of the integrand in 
Eq. (jl9p . allows to simultaneously compute the moments 
dAdijk (N, M) for all N, M, without having to do an 
individual Monte-Carlo integration for each of them as 
in ([TO]) . This greatly reduces the computational bur- 
den. To give an example, for our global analysis of polar- 
ized PDFs, which currently includes about 50 very time- 
consuming NLO calculations of single-inclusive jet and 
pion production at RHIC in each step of the mini- 
mization, all the necessary grids in Eq. (j23p can be com- 
puted within approximately one day of CPU time on a 
single standard PC with a CoreDuo processor running 
at 2 GHz. Once these grids are available, a full PDF 
fit can be finalized in about 15-30 minutes. A detailed 
PDF uncertainty assessment, which requires a very large 
number of minimizations, can then be performed eas- 
ily in about 1-2 weeks. We note that in practice one does 
not even need to physically store the set x^j^\ X2^\ w^J^', 
since the summations (1231) can in fact be done simulta- 
neously with the generation of the events n. 

The formula in Eq. (|20[) applies to any computation 
of the theoretical cross section when a Monte-Carlo inte- 
gration is employed. This is the case for analytical NLO 
calculations of single- inclusive jet or hadron rates [53l.[53| 
where Monte-Carlo integration techniques are just used 
to perform the convolution with the parton distributions, 
or for NLO parton-level Monte-Carlo generators evaluat- 
ing general infrared-safe observables [55, 56] . Therefore 
our method based on Mellin moments and Monte-Carlo 
sampling techniques outlined in Eqs. (fT51) - (P51) above is 
completely general and can be straightforwardly applied 
to any observable for which a perturbative QCD descrip- 
tion is available. It can be applied to global analyses 
of polarized PDFs, pursued in this pap er, extractions of 
fragmentation functions, see Refs. [37|, but equally well 
to analyses of ordinary spin-averaged PDFs incorporat- 
ing Tevatron and future LHC data consistently at NLO 
or beyond. We note that for the latter case an alterna- 
tive method, called "fastNLO" 57], has been developed, 
which also allows to include NLO corrections to hadronic 
scattering in a fast and efficient way in a global PDF 
analysis. Like our method, it amounts to preparing huge 
look-up tables that contain all the time-consuming NLO 
calculations prior to the actual fit. In the "fastNLO" 
case, this is done in a;-space, and interpolations between 
the various support points in x are used. This step ap- 
pears to be rather time-consuming and the method has 
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so far been tested only for inclusive jet production data. 
In addition, the evolution of the PDFs needs to be per- 
formed as a separate calculation, whereas in our approach 
it is immediately included in Mellin moment space. 

This latter advantage can in fact be used for further 
improvements of our Mellin technique. The scale depen- 
dence of the PDFs, which we have so far suppressed, can 
be schematically written as 

AMN,^i) = Y,E^HN,^lo,^^)AUiN,^lo). (24) 

i' 

Here, Eiii{N, /iq, fj,) denotes the appropriate evolution 
matrix from the initial scale /iq where we parameterize 
the PDFs to the scale relevant for a certain observ- 
able O. Inserting into Eq. (fT7|) allows to absorb also 
the scale evolution of the PDFs into the pre-calculated 
Mellin grids by extending Eq. ([23]) to 

dA^,,fe(iV,M) = -VVa:r*'x2-^ 

n—1 i'j' 

Ee,{N, fIo,^l)E,,,{M, a^o, ^ ,(25) 

so that the luminosity function in Eq. psp now only 
refers to the PDFs at the initial scale /io, L.ij{N,M) — 
Afi{N, /ip) A/j(Af, fio). An advantage of this re-shuffling 
is that now all dependence on the scale fi is contained 
in the look-up tables eliminating the need to per- 

form the scale evolution later in the fitting code. More 
importantly, if the experimental observable used in the 
global fit involves an integration over a bin in, say, the 
transverse momentum px of a jet, which also acts as the 
factorization scale, the n dependence of the PDFs is cor- 
rectly taken care of in the integration. While we have not 
made use of this particular improvement in the present 
analysis, we expect it to be useful in the future for further 
optimizing the performance of the global analysis code. 

III. GLOBAL ANALYSIS AND UNCERTAINTY 
ESTIMATES FOR POLARIZED PDFS 

In this Section we give a detailed account of the first 
global analysis of polarized PDFs presented in Ref. 
which in the following will be referred to as "DSSV" . We 
first discuss the data selection and the determination of 
the best fit, which we compare to the fitted data. We then 
focus on the studies of uncertainties, including a compar- 
ison of the Lagrange multiplier method used in [ll] and 
the more standard Hessian error matrix approach. For 
the latter we present a new family of eigenvector PDFs, 
as described above, which greatly facilitates estimates of 
the PDF uncertainties of any given observable of interest. 

A. Determination of the optimal fit 

Our first physics objective is to establish the set of 
polarized PDFs that gives the optimum theoretical de- 



TABLE I: Data used in our NLO global analysis of polarized 
parton densities, the individual values for each set, and 
the total of the fit. 



experiment 


process 


Adata 




EMC [2] 


BIS (p) 


10 


3.9 


SMC [3] 


BIS (p) 


12 


3.4 


SMC 13} 


BIS (d) 


12 


18.4 


COMPASS [4] 


BIS (d) 


15 


8.1 


E142 [5] 


BIS (n) 


8 


5.6 


E143 [6] 


BIS (p) 


28 


19.3 


E143 \6} 


BIS (d) 


28 


40.8 


E154 [7] 


BIS (n) 


11 


4.5 


E155 ]8} 


BIS (p) 


24 


22.6 


E155 [9] 


BIS (d) 


24 


17.1 


HERMES [10] 


BIS (He) 


9 


6.3 


HERMES [llj 


BIS (p) 


15 


10.5 


HERMES [11] 


BIS (d) 


15 


16.9 


HALL-A [12] 


BIS (n) 


3 


0.2 


CLAS [13| 


BIS (p) 


10 


5.9 


CLAS [13] 


BIS (d) 


10 


2.5 


SMC [14] 


SIBIS {p,h+) 


12 


18.7 


SMC [14] 


SIBIS {p,h-) 


12 


10.6 


SMC [14| 


SIBIS {d,h+) 


12 


7.3 


SMC [14] 


SIBIS {d,h-) 


12 


14.1 


HERMES [15] 


SIBIS (p,ft+) 


9 


6.4 


HERMES [15] 


SIBIS (p,/i~) 


9 


4.9 


HERMES [15] 


SIBIS (d,/i+) 


9 


11.4 


HERMES [15] 


SIBIS (d,/i-) 


9 


4.5 


HERMES [10] 


SIBIS (He,/i+) 


9 


4.7 


HERMES [10] 


SIBIS {Re,h ) 


9 


6.9 


HERMES [15] 


SIBIS (p,7r+) 


9 


9.6 


HERMES [15] 


SIBIS (p,7r-) 


9 


4.9 


HERMES [15] 


SIBIS (d,7r+) 


9 


9.4 


HERMES [15] 


SIBIS (d,7r-) 


9 


19.5 


HERMES [15] 


SIBIS {d,K+) 


9 


6.2 


HERMES fl5] 


SIBIS {d,K-) 


9 


5.8 


HERMES [15] 


SIBIS {d,K++K-) 


9 


3.4 


COMPASS [16] 


SIBIS (d,/i+) 


12 


6.2 


COMPASS [16] 


SIBIS (d,/i-) 


12 


12.0 


PHENIX [22] 


pp (200GeV,7r") 


10 


14.2 


PHENIX [23] 


pp (200 GeV, tt") 


10 


7.1 [13.8] " 


PHENIX [24] 


pp (62 GeV, TT°) 


5 


3.1 [2.8] " 


STAR [25] 


pp (200 GeV, jet) 


10 


8.8 


STAR (prel.) [26] 


pp (200 GeV, jet) 


9 


6.9 


TOTAL: 




467 


392.6 



"The PHENIX data were still preliminary when the global anal- 
ysis presented here was performed. The value quoted in 
brackets is evaluated with the published data |23| . |24| but without 
re-fitting. 



scription of the available hard scattering data. The data 
sets for the spin asymmetries we use in our analysis are 
listed in Tab. HI along with the number of data points 
fitted. We minimize the effective function in Eq. ([1]). 
Attempts to further improve the global fit by introducing 
normalization shifts for each experiment and minimizing 
X^ according to Eq. ^ were to no avail. All theoretical 
spin asymmetries in Eq. ^ are calculated at NLO, using 
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the appropriate factorized leading- twist expressions. We 
use the MS scheme throughout, and all our results for 
the polarized PDFs will refer to this scheme. 

In case of inclusive DIS, the asymmetries are com- 
puted, as in our previous analyses [3l|,[3^, as the ratios 
between the polarized and unpolarized structure func- 
tions, 

with 

+ ^[AC,®{Aq + Aq) + ACg®Ag]} , (27) 

and the corresponding expression for Fi{x,Q^), both 
computed at NLO using the appropriate coefficient func- 
tions [13, [111. For DIS off a deuteron target, we take into 
account the ujd = 5.8% D-wave state probability in re- 
lating the gi structure function of the deuteron to those 
of proton and neutron: gf = {1 — 1.5ujd){9i+9i)/'2'- The 
extension of the expression in ([^5]) to SIDIS is straightfor- 
ward, using the NLO coefficient functions given in [59| . 

For the case of pp scattering, the spin asymmetries 
are computed using the generic expression in Eq. pH]) at 
NLO, and its spin- averaged counterpart. The NLO cor- 
rections for high- pt single- jet and hadron production can 
be found in [53l |5J] , respectively. We have always cho- 
sen the renormalization and factorization scales as the 
transverse momentum pT of the observed final state, a 
choice that leads to good agreement of NLO calculations 
[13, [H, m and experimental data from RHIC [H |6^ 
in the spin-averaged case. For the computation of the 
unpolarized cross sections, we always use the NLO unpo- 
larized PDFs of Ref. [131 • Whenever fragmentation func- 
tions are needed, as is the case for SIDIS and the RHIC 
pp nX data, we use the DSS set ^37*1 for pions, kaons, 
and unidentified charged hadrons, which was recently ob- 
tained from a global analysis of hadron production data. 
The use of up-to-date fragmentation functions that are 
consistent with HERMES and RHIC unpolarized data 
and have quantified uncertainty estimates [S^l is a cru- 
cial ingredient of our analysis. It is a major difference 
with respect to the polarized PDF analysis in Ref. [36j . 
where also SIDIS data were included. In the computa- 
tion of the contribution from SIDIS asymmetries we 
have taken into account the uncertainty coming from the 
set of fragmentation functions. In practice, this was done 
by determining for each data point the maximum varia- 
tion of the theoretical estimates Tj in Eq. ([T]) due to the 
FFs within their own uncertainty ranges quoted in [37l|. 
This variation was treated as an additional uncertainty 
and added in quadrature to the experimental error SDj. 
Since at least for pions the FFs are fairly well determined, 
this amounts to additional uncertainties of a few percent 
at most. We do note, however, that this uncertainty es- 
timate does not reflect possible systematic problems in 



interpreting the available spin-averaged kaon SIDIS data 
within a leading-twist pQCD analysis, a point to which 
we will return later. 

We do not systematically include any higher twist (or, 
more generally, power-suppressed) contributions in the 
theoretical calculations that enter our analysis, neither 
for the inclusive and semi-inclusive DIS observables, nor 
for proton-proton collisions. We employ a cut of > 
1 GeV^ for all DIS and SIDIS data, and of pr > 1 GeV 
for the RHIC high-p^ polarized pp data. These cuts have 
the purpose to exclude regions where contributions be- 
yond the leading twist, factorized framework of pQCD 
become crucially important. For example, for the SIDIS 
and RHIC pp data, it is known that the underlying un- 
polarized cross sections in the same kinematic domain, 
i.e., for scales above 1 GeV, can be quite successfully 
described by pQCD [sj. 

That said, it is known [6l| that the relation between Ai 
and gi/Fi in Eq. (^5]) is corrected by a factor (1-1-7^) = 
(1 -I- 4M^x^/Q^) on the right-hand-side (r.h.s.), corre- 
sponding to a target mass correction. It has been pointed 
out [6^ that this correction is non-negligible in some 
kinematic regimes accessed by the lower energy fixed- 
target experiments, typically at relatively low and 
high X. We have therefore corrected for this factor where 
necessary. Specifically, since our set of polarized PDFs is 
defined and related to the measured asymmetries through 
the leading- twist relations ([26]) and ((27)) . our choice is 
to multiply data sets that are given in terms of mea- 
sured gi/Fi by the factor (1 -I- 7^), but to leave data 
sets for measured Ai unchanged. The resulting data 
are confronted with the NLO leading twist calculation. 
We stress that various other choices have been adopted 
in the literature [H, [s^ [13], using, for example, pa- 
rameterizations of experimental data for F2{x,Q^) and 
R{x,Q^) = F2/{2xFi) — 1. All choices are equivalent in 
the strict leading twist sense, but will in general differ in 
the amount of power corrections needed to describe the 
data at lower and/or higher x. As we shall see be- 
low, we find that for our choice without inclusion of any 
power corrections an excellent description of all data sets 
within our specified cuts is achieved (cf. also Tab.|T|), with 
no visible discrepancies even at rather low scales. For the 
time being, we thus regard our choice as one that "em- 
pirically" alleviates the need for power-suppressed cor- 
rections in spin asymmetries. Despite some significant 
progress on the analysis of higher twist effects in polar- 
ized DIS [H, [63 , further detailed investigations will be 
needed in this area, including implementation of the full 
target-mass corrections [i^, the effects of yet higher or- 
ders [13, [sBl which generally reduce the need for higher 
twist contributions [M, iSQ], and so forth. 

In order to find the optimal helicity-dependent PDFs 
from a minimization, we parameterize them at an in- 
put scale of /Lto = 1 GeV by the following flexible form [28j 

xAMx,fil)^N,x'''il~xf^{l+j^V^ + Tj,x), (28) 

with independent parameters for Au + Au, Ad + Ad, 
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Au, Ad, As = As, and Ag (note that here and in the 
following we interchangeably use A/„ = Au, Afg — Ag 
etc. to denote the polarized PDFs). The minimiza- 
tion is carried out with respect to the set of parameters 
{at} = {Ni,ai, f}i,^i,rii}. The PDFs are evolved to the 
scales /i > /io relevant in experiment. The particular 
functional form and the value for fiQ are not too cru- 
cial, as long as the parameterization is flexible enough 
to accommodate all hard scattering data within their 
ranges of uncertainties. The ansatz ipSj) deviates consid- 
erabl y fr om the standard form used in fits to DIS data 
only |3ll . [3^ . [33I [13, [11], where = r]i = 0, inasmuch as 
it allows the PDFs to develop nodes and to deviate from 
an SU(3) flavor symmetric sea. As will be seen from our 
results presented below, this extra freedom in parameter 
space {oi} is crucial in a comprehensive analysis of DIS, 
SIDIS, and RHIC pp data. 

In addition to the much more flexible input parame- 
terization proposed in the preceding paragraph, we have 
repeated the analysis with alternative parameterizations, 
some of them even more flexible than the one we choose. 
For example, we have chosen the powers of x in the last 
two terms different from 1/2 and 1 , even allowing the fit 
to vary them. None of these modifications resulted in 
any significant improvement in the quality of the fit to 
data, or changes of the uncertainty bands. This indicates 
that the present data is not really able to discriminate 
between various forms of the input distributions, as long 
as a sufficiently flexible choice is made. Therefore our 
present choice does not introduce large additional uncer- 
tainty in that respect. 

Analyses of polarized PDFs routinely use constraints 
that can be derived from baryonic semi-leptonic /3-decays 
under the assumption of SU(2) and SU(3) flavor symme- 
tries [S^l • These relate combinations of the first moments 
of the PDFs to the constants F and D parameterizing the 
/?-decay rates. We make use of these constraints in our 
present analysis; however, rather than imposing the ex- 
act SU(2) and SU(3) flavor symmetry relations, we allow 
for deviations in our fit within the uncertainty ranges of 
the F,D values. Specifically, we set 

AE„-AI]<i - (F + i?)[l + esu(2)], (29) 
AI]„ + ASd~2AI], = (3F-i?)[l + £su(3)], (30) 

where 

AS/ ^ [Afl + Afl] (m^) ^ [A/, + A/,] {x, t^l) dx , 

Jo 

(31) 

esu(2,3) parameterize the departures from exact SU(2) 
and SU(3) symmetries, and where we use the latest values 
F + D = 1.269 ± 0.003 and 3F - D = 0.586 ± 0.031 (see, 
e.g., Ref. (Tlj). As a practical matter, we trade the input 
parameters N^+u and N^j^^ in Eq. (|28p for esu(2,3) and 
fit the latter. Here the relative uncertainties of -I- D 
and 3F — D are assumed to represent the typical ranges 
of esu(2,3); we use them to include the deviations from 
esu(2,3) = as additional contributions to 1 similarly to 



the case of normalization uncertainties shown in the first 
term of Eq. ^ . We note that the relative uncertainties 
oi F + D and 3F~D are rather modest and may not fully 
refiect the actual breaking of the SU(2) and, in particular, 
SU(3) symmetries, for which larger breaking effects have 
been discussed in the literature [68| . This issue may need 
to be revisited in the future. For now we note that as a 
result of this the PDFs in our fits will naturally have a 
tendency to have relatively small esu(2,3)- 

Rather than determining also the strong coupling as 
in the global fit along with the PDFs, we take Aqcd = 
334.2 MeV for n/ = 4 flavors from the analysis of unpo- 
larized PDFs in Ref. [43|. The scale dependence of as 
is computed by numerically solving its renormalization 
group equation at NLO accuracy. The charm and bot- 
tom quark thresholds are crossed at nic = 1.43 GeV and 
TOfc = 4.3 GeV, respectively. As already mentioned, the 
scale evolution equations for the PDFs are solved analyti- 
cally in Mellin moment space by explicitly truncating the 
solutions at NLO. Likewise, all observables used in our 
fit are computed consistently at NLO accuracy in the MS 
factorization scheme. All quarks are treated as massless, 
and charm and bottom PDFs are turned on in the evolu- 
tion at Q = mc,b- We note that for all presently available 
spin-dependent observables heavy quarks play a negligi- 
ble role, and, for the time being, we can safely refrain 
from introducing more elaborate variable flavor number 
schemes which take into account quark mass effects near 
threshold (see, e.g., Refs. [29l.[30|). 



TABLE II: Parameters {a^} describing our optimum NLO 
(MS) xAfi in Eq. ([28} at the input scale ^io = 1 GeV. 



flavor i 




Oii 




7i 




u + u 


0.677 


0.692 


3.34 


-2.18 


15.87 


d+d 


-0.015 


0.164 


3.89 


22.40 


98.94 


u 


0.295 


0.692 


10.0 





-8.42 


d 


-0.012 


0.164 


10.0 





98.94 


s 


-0.025 


0.164 


10.0 





-29.52 


g 


-131.7 


2.412 


10.0 





-4.07 



The parameters {a^} representing our best global fit of 
polarized parton densities A/j in Eq. (p8|) , henceforth de- 
noted as "set iS*"" , are given in Tab. |TT1 A few additional 
remarks are in order here. The currently available data 
do not fully constrain the entire x dependence of Afi 
imposed in Eq. (|28p . and we are forced to make some re- 
strictions on the parameter space {a,}. For the sea quark 
and gluon densities we set 7^ = in Eq. ([28]) . This only 
marginally limits the freedom in the functional form and 
still allows nodes. In addition, we tie the small x behav- 
ior, represented by the ai in Eq. (pS)) . of Au -f Au and 
Ad -I- Ad to that of Au and Ad, respectively, which is 
reasonable as sea quarks likely dominate in this region. 
The parameter ag always came out close to a^, so we 
set them equal. Since the SIDIS data are not yet suffi- 
cient to distinguish As from As, we assume As = As 
throughout. In general, the x ^ 1 behavior of all PDFs 
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is rather unconstrained as is the case even for current 
sets of unpolarized PDFs. This is because there are no 
data sensitive to a; > 0.6. To avoid problems with the 
fundamental positivity constraint 

|dAcr| < da , (32) 

we make sure that all A/^ vanish at large x at least as fast 
as the fi of our reference set of unpolarized PDFs 
which puts constraints on the (3i. Choosing any other 
recent set of unpolarized PDFs, like CTEQ6 [4^, does 
not alter our results. We note that implementing the 
positivity constraint at the level of PDFs [6§|, i.e., 

\AMx,Q^)\<Mx,Q') (33) 

is strictly valid in LO only, but in the MS scheme is 
sufficient to guarantee ([5^ also in NLO. The parame- 
ters Pu+u and come out very close to their cor- 
responding values for the unpolarized case, implying 
{Au -t- Au)/{u + u) ^ const, as a; — s- 1, and likewise 
for (Ad + Ad)/{d + d). Since the other (3i in Eq. ([28]) 
are only very weakly determined by the fit, we fix them 
to their preferred values within the positivity constraints 
whenever we examine uncertainties of the PDFs, in or- 
der to avoid extremely flat directions in parameter space 
where varies only very slowly. Therefore, and as is 
in general the case in PDF studies, our uncertainty esti- 
mates for helicity-dependent PDFs are valid only in the 
x-region explored by experiment. Notice that the near 
equality of f?(j+ j and Tyj is not imposed but a result of the 
fit (in fact, the actual values for these parameters before 
rounding are 98.9384 and 98.9354, respectively). 

In total this leaves us with 19 free parameters in the 
fit [or 17, if we fix also pu+u and which we include 

later on also in our uncertainty estimates. We tried to 
relax the imposed constraints discussed above, but found 
that present data, i.e., the effective function, are not 
really sensitive to them. In Tab.HIlwe have converted the 
fitted values for £su(2,3): defined in Eqs. ((^ and 
back to Nu+u and A^^+j for convenience. For the optimal 
DSSV fit (set 5°) we find 

£su(2) = 0.0011 , esu(3) = -0.0035 , (34) 

which corresponds to only very minor violations of the 
canonical constraints on the first moments AE„ — AE^ 
and AT,u + AE^ — 2AEs assumed in most fits so far. As 
we have discussed above, the smallness of esu(2,3) is not 
really a surprise in view of the relatively small nominal 
uncertainty of the F + D and 3F — D values in Eqs. ([^5]) 
and ((30)) . If correct, the small value for esu(3) has inter- 
esting implications on the behavior of the As(a;) = As{x) 
distribution in the best fit, as we shall see later. 

B. Comparison to fitted data 

The total of the best fit S" is 392.6 for 467 data 
points used in our NLO global analysis. We list in Tab.|T] 



also the individual values for each experiment. As one 
can see, there are only very few cases where the X^/^d^ta 
is significantly larger than one. In each case, this is due 
to large fluctuations of some of the data points in that 
particular set which are impossible to accommodate in 
the fit. Figure [T] shows the comparison of our fit to the 
fitted DIS data, while the comparison to the SIDIS data 
is shown in Fig. [2l Notice that in Fig. [T] the plots are 
generically labeled as asymmetries "^i"; however, in the 
case of the E143, E155, CLAS and Hall A data, they ac- 
tually show the reported structure function ratios and are 
compared to the DSSV estimates for the asymmetries, di- 
vided by the factor (l-|-7^). The overall agreement of the 
experimental sets in the global analysis is excellent. All 
data can be very satisfactorily described by a universal 
set of polarized PDFs as is assumed by the fundamental 
factorization theorem. The agreement with the RHIC pp 
data is equally good; we have shown the corresponding 
comparison in our previous paper [28i] and will come back 
to it in the next Subsection. 

Figures [T] and [2] also show the results obtained for the 
set of polarized PDFs of 36] in the following labelled as 
"DNS". Apart from the fact that not all of the present 
data sets were available at the time of ,36 j , a main differ- 
ence between the two analyses resides in the fragmenta- 
tion functions used when including the SIDIS spin asym- 
metries in the fit. To illustrate this point, we have used 
here the new fragmentation functions of [33] also for the 
calculations with the DNS set [36]. As can be seen from 
Fig. [21 this leads to significant differences, in particular, 
in the kaon and, to a lesser extent, the isospin sensitive 
(proton target) SIDIS asymmetries. This is mainly due to 
the strange fragmentation functions, which directly affect 
the strange quark polarization, and also to differences in 
the light sea quark distributions. Figure [1] shows that 
there is, however, little difference between the two sets 
as far as the inclusive DIS asymmetries are concerned. 
The changes in the strange quark and other sea polariza- 
tions are thus compensated here, either mutually or by 
the other parton distributions. 



C. Extracted PDFs, their uncertainties, and their 
physics 

Figure [31 shows the extracted polarized PDFs at = 
10 GeV^, along with estimates of their uncertainties for 
the Hessian and Lagrange multiplier methods, both for 
a tolerance of Ax^ = 1. The results for the Lagrange 
multiplier method were already shown in our previous 
paper [2^ along with a more conservative estimate of 
the PDF uncertainties based on Ay^ jy^ — 2%. For this 
method, the estimates were obtained by varying the first 
moments of the distributions, truncated to the region of 
momentum fractions 0.001 < x < 1 covered by the data 
included in the fit. The distributions corresponding to 
the maximum variations of the truncated moments for a 
given increase Ax^ were then taken as faithful estimates 
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FIG. 1: Inclusive DIS spin asymmetries 0, S i, i, i, Q, i, i, E, 111111' El compared to the best fit results of our global 
analysis ("DSSV", solid lines), and for the set of polarized PDFs of M)\ ("DNS", dashed lines). 



of the range of variation of the PDFs. In the case of the 
polarized gluon distribution, this procedure was found 
to be not adequate [1^ because of the fact that there 
is a significant amount of rather precise proton-proton 
colhsion data constraining the gluon density in how- 
ever a relatively narrow region of momentum fraction, 
0.05 < a; < 0.2. In this way the variations of the gluon 
distribution's integral in the full region 0.001 < a; < 1 
tend to produce distributions that favor the variations 
outside the pp kinematic region, misrepresenting the un- 
certainties inside. In order to circumvent this problem, 
we performed variations of the integral of the gluon dis- 
tribution in three different x regions, 0.001 < a; < 0.05, 
0.05 < X < 0.2, and 0.2 < x < 1, allowing them to jointly 
contribute a change in of Ax^ = 1 . Clearly, the choice 
for these regions and the way they share the increase in 
Ax^ is not unique. In order to specifically focus on the x- 
region accessed at RHIC, we also performed a dedicated 
study of the truncated moment of in this region, al- 
lowing variations of Ax^ = 1 from this region alone. The 
results for the truncated moments of our polarized PDFs, 



TABLE III: Truncated first moments A/j'['' °°^^^' at = 
10 GeV^ and their uncertainties for Ax^ ~ 1 obtained with 
the Lagrange multiplier and the Hessian methods. For future 
reference, we also recall the results for the Lagrange multiplier 
method obtained in [28l ] under the assumption Ax^/x^ = 2%, 
which are to be considered more realistic estimates of the 
uncertainties. In the last line, Ag^™'^ represents the first 
moment but truncated to [0.05 0.2]. 

Lagr. Ax^ = 1 Hessian Lagr. Ax^/x^ = 2% 



Au + Au 


n 7QQ+0011 


0.793±0.012 


Q 7qQ+0.028 
U. ( »0_g 034 


Ad + Ad 


-0.416«;S^ 


-0.416±0.011 




Au 


0.028tH^J 


0.028±0.022 


0.028tHf9 


Ad 


-0.089««^^ 


-0.089±0.029 




As 


-0.006«:«1» 


-0.006±0.012 


-o.ooelH^? 


AS 


U.aDD_Q 018 


0.366±0.017 


U.JDD_o.o62 


Aff 


U.UiO_o 120 


0.013±0.182 


o.oisi":™^ 




n oots+oo^i 

U.UUO_o Q5g 


0.005±0.056 


o.oo5l«:lf. 



AMx,Q^)dx, (35) 



are given in Tab. IIIII 

Inspection of Fig. [3] and Tab. IIIII reveals that the Hes- 
sian and the Lagrange multiplier methods yield fairly 
similar Ax^ = 1 uncertainties, except for the spin- 
dependent gluon distribution, for which the Lagrange 
multiplier approach yields a still significant but gener- 
ally smaller uncertainty than the one predicted by the 
Hessian method using Eq. PT|) . The agreement between 
the two often becomes better when the observable is bet- 
ter constrained by the data, as is the case for the integral 



of Ag over only the x-range probed at RHIC, or for the 
actual physical observables that determine A^. As an 
example, in Fig. 3] we show the estimated uncertainties 
for the double-longitudinal spin asymmetry, 



(7 ' ' — (7 



(36) 



for pp tt'^X at RHIC, where the superscripts denote 
the helicities of the incoming protons, computed with 
both the Lagrange multiplier and the improved Hessian 
approaches. As can be seen, the two give very similar re- 
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FIG. 2: Same as Fig. [T] but for the semi-inclusive DIS asymmetries [lol . [T3 . fisl . [l6l ]. In all calculations the fragmentation 
functions of have been used. 



suits. This feature can be traced back to correlations be- 
tween the parameters, in the sense that some of them can 
compensate variations forced in the others. We note that 
such kinds of correlations are fully accounted for in the 
Lagrange multiplier approach, whereas it is not generally 
clear how well are they represented by the approximated 
Hessian matrix. We shall investigate the distinctive fea- 
tures between the two methods later, but will focus first 
on the physics aspects related to our extracted polarized 
PDFs. 

Table IIVI shows the evolution of the central values for 
the truncated first moments /^f^'^^-^^^^^^ with Q^. AS 
denotes the quark singlet combination, i.e., the sum of all 
quarks and anti-quarks. We also show the evolution of 
the full first moments Af^. These obviously rely on an 
extrapolation of the PDFs to x-values outside the mea- 
sured region, and it is difficult to estimate the uncertainty 
associated with this. 

Total up and down distributions: Au + AU and Ad -t- Ac?, 
which inclusive DIS probes primarily, are the by far best 
determined distributions. Their uncertainty bands are 



very narrow, see Fig. [31 and also our results agree very 



32, 
70| 



well with the determinations in previous analyses |3ll . 
SHiil. We note that recent lattice QCD results 
of the full first moments AI]„ = Au^ + Au^ and AE^ = 
Ad^ + Ad^ (albeit excluding disconnected diagrams) also 
agree very well with the values we extract, which may 
shed light on the validity of assumed extrapolations of 
the parton distribution functions to small x. 

We have mentioned earlier that in our fit i?„ = (Am -|- 
Au)/{u + u) and Rd = {Ad + Ad)/{d + d) become con- 
stant in the "valence region" as x — > 1, where the sea 
quark contributions become small. Figure [5] shows the 
ratios Ru, Rd along with the most relevant experimental 
data. The information at the highest values of x comes 
from the Jefferson Laboratory Hall- A experiment p^. 
As one can see, our i?„ goes to unity at high x, which 
is consistent with expectations in relativistic constituent 
quark models [71[, but also in perturbative QCD, using 
power counting and hadron helicity conservation |72l | . We 
furthermore find that Rd remains negative in the region 
where it is constrained by data and presently shows no 
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TABLE IV: Truncated first moments, A//'[° '"'^"*^' , and full ones, Afl, of our polarized PDFs at various Q^. 



x-range in Eq. (I35II 


[GeV^] 


Au + Au 


Ad + Ad 


Am 


Ad 


As 


As 


AE 


0.001-1.0 


1 


0.809 


-0.417 


0.034 


-0.089 


-0.006 


-0.118 


0.381 




4 


0.798 


-0.417 


0.030 


-0.090 


-0.006 


-0.035 


0.369 




10 


0.793 


-0.416 


0.028 


-0.089 


-0.006 


0.013 


0.366 




100 


0.785 


-0.412 


0.026 


-0.088 


-0.005 


0.117 


0.363 


0.0-1.0 


1 


0.817 


-0.453 


0.037 


-0.112 


-0.055 


-0.118 


0.255 
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0.814 


-0.456 


0.036 


-0.114 


-0.056 


-0.096 


0.245 




10 
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-0.458 


0.036 
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-0.057 


-0.084 


0.242 




100 


0.812 


-0.459 
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-0.116 


-0.058 


-0.058 
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-0.04 
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FIG. 3: Our polarized PDFs of the proton at = 10 GeV^ 
in the MS scheme, along with their Ax^ = 1 uncertainty 
bands computed with Lagrange multipliers and the improved 
Hessian approach, as described in the text. 



tendency to turn towards -1-1 at high x. The latter be- 
havior would be expected for the pQCD based models. 
We note that it has recently been argued [t^ that the 
upturn of Rci in such models could set in only at rela- 
tively high X, due to the presence of valence Fock states of 
the nucleon with nonzero orbital angular momentum that 
produce double-logarithmic contributions ^ ln^(l — x) in 
the limit of x —^ 1 on top of the nominal power behav- 
ior. The corresponding expectation is also shown in the 
figure. In contrast to this, relativistic constituent quark 
models predict Rd to tend to —1/3 as a: —> 1, perfectly 
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-0.002 



-0.004 




m!i$ Ax - 
Ax': 



:1 (Lagr. multiplier) 
:1 (Hessian) 



FIG. 4: Uncertainties of the calculated at RHIC in our 
global fit, computed using both the Lagrange multiplier and 
the Hessian matrix techniques. We also show the correspond- 
ing PHENIX data [13]. 

consistent with the present data. 

Light sea quark polarizations: The light sea quark and 
anti-quark distributions turn out to be better constrained 
now than in previous analyses 
of more precise SIDIS data 0,1 



36 I , thanks to the advent 
14, 111, Ell and of the new 
set of fragmentation functions [37[ that describes the ob- 
servables well in the unpolarized case. Figure [H] shows the 
changes in of the fit as functions of the truncated first 
moments Au^-^°-°°^^^\ Ad^'^°-°°^^^^ defined in Eq. ([35|l. 
obtained for the Lagrange multiplier method. On the 
left-hand-side, Figs. [S] (a), (c), we show the effect on the 
total x^, as well as on the values for the individual 
contributions from DIS, SIDIS, and RHIC pp data and 
from the F, D values. It is evident that the SIDIS data 
completely dominate the changes in x^- On the r.h.s. of 
the plot, Figs.[n](b), (d), we further split up Ax^ from 
SIDIS into contributions associated with the spin asym- 
metries in charged pion, kaon, and unidentified hadron 
production. One can see that the latter dominate, closely 
followed by the pions. The kaons have negligible impact 
here. For Au^'t" °"^^^l , charged hadrons and pions are 
very consistent, as far as the location of the minimum 
in x^ is concerned. For Aj^'t^ *^"^^^! there is some slight 
tension between them, although it is within the tolerance 
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FIG. 5: High- a; behavior of our Ru = (Ait + Au)/{u + u) 
and Rd = (Ad + Ad)/{d + d), shown by the solid curves, 
along with the data from p^. [T5| . as shown in [t^I- The 
dashed lines present the predictions based on power counting 
and perturbative QCD, taking into account nonzero orbital 
angular momentum Fock states 73j. We note that these re- 
sults, as well as some of the experimental data, are for the 
ratios Au/u and Ad/d rather than for Ru and Rd, which, 
however, makes little difference for the large x values we con- 
sider here. The arrows show the expectations for Au/u and 
Ad/d in relativistic constituent quark models [tH ]. 




^^1, [0.001^1] ^^1, [0.001^1] 

FIG. 6: (a) The profiles and partial contributions Ax^ of 
the various types of data sets for variations of the truncated 
first moment Au^-^°-°°^^^-°^ at = 10 GeV^ (b) Partial 
contributions Ax^ for the various hadrons produced in SIDIS. 
(c) and (d): same as (a) and (b), but for Ad^-^°-°°^^^-°l 



of the fit. 



Of particular physics interest is a possible flavor sym- 
metry breaking in the light sea, i.e., Au ^ Ad, given the 
well-established significant difference between u and d in 
the spin-averaged case @, [s^l- Figure [3] indeed clearly 
points to a largely positive Am distribution, but a neg- 
ative (and larger) Ad. Figure [7| specifically shows the 
difference x{Au — Ad), which is positive within uncer- 
tainties. Note that we show both the Ax^ — 1 and the 
more conservative Ax^/x^ = 2% uncertainty bands here. 

The pattern of symmetry breaking in the light anti- 
quark sea polarizations shown by Figs. [3] and [7] has been 
predicted at least qualitatively by a number of models 
of nuclcon structure. A simple intuitive consideration of 
the Pauli principle roughly gives the observed picture: 
if valence-u quarks primarily spin along the proton spin 
direction, uu pairs in the sea will tend to have the u 
quark polarized opposite to the proton. Hence, if such 
pairs are in a spin singlet, one expects Am > and, by 
the same reasoning, Ad < 0. Expectations based on the 
Pauli principle have been made quantitative in [tII and 
the "valence" scenario of [3l|, and the resulting predic- 
tions are shown by the dot-dashed line in Fig. [71 They 
tend to lie somewhat higher than our extracted Am — Ad, 
but are certainly qualitatively consistent, given the still 
relatively large uncertainties. The same is true for the 
case of the chiral quark soliton model W^, represented 
by the dotted line in the figure. Within the large-iVc 
limit of QCD on which this model is based, one in fact 
expects I Am — Ad| > |m — d|. As comparison of our ex- 
tracted x{Au — Ad) with the result of [46] for x{d — u) in 
Fig. [7] shows, one can presently not yet decide whether 
this expectation is fulfilled. Predictions for Am — Ad 
have also been obtained within meson cloud models [t^ ; 
it has been found in [T^ that also here a flavor asym- 
metry of similar size is possible. Finally, also statistical 
parton models [H, ItI] obtain a similar size of Am — Ad. 
We note that predictions for the individual Am and Ad, 
where available, agree on Am > 0, Ad < 0, consistent 
with our results in Fig. [31 but may differ in the rela- 
tive size of the distributions. For example, the results 
of [m, [zil have |Ad| > Am, as in Fig. [3 while the sta- 
tistical models find the two distributions to be of more 
equal absolute size. 

Strange quark polarization: The polarization of strange 
quarks has been a focus since the very beginning of the 
proton spin crisis. The reason is that in the parton model 
and assuming SU(3) symmetry (see Sec. IIII A[) one has 

AS = E„ + Ed + = (3F -D) + SAS^ , (37) 

where the AE/ are as defined in Eq. (|3T|) but now for 
arbitrary scale Q, and AE is the total quark and anti- 
quark spin contribution to the proton spin. If the latter 
is found to be small experimentally, AE ~ 0.25, the im- 
plication is that strange quarks make a significant nega- 
tive contribution to the proton spin. Indeed, most fits to 
only inclusive DIS data have preferred a large and nega- 
tive strange quark polarization. The same was found in 
Ref. [3^, even though here the SU(3) flavor symmetry 
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FIG. 7: The difFerence between xAu and xAd at = 
10 GeV'^, along with the uncertainty bands for Ax^ ~ 1 and 
Ax^/x^ — 2%. The dot-dashed and dotted lines show the pre- 
dictions of the valence scenario of [3l| ] and the chiral quark 
soliton model of [t^ , respectively. We also show the result ob- 
tained in an earlier global analysis fsH of DIS and SIDIS data 
(light dotted), for which the fragmentation functions of ^3] 
were not yet available. The dashed line displays for compari- 
son the flavor asymmetry x{d — u) in the spin-averaged case, 
using the PDFs of 



was not enforced. 

At variance with these results, the best fit in our 
present analysis has a polarized strange distribution As 
that is positive at large a;, but negative at small momen- 
tum fractions. Before we discuss the origin and signifi- 
cance of this result, we note that a prerequisite for it is 
that we have adopted a more flexible parameterization for 
the strange quark distribution in this work, which per- 
mits a node. This is again in contrast with the previous 
fits in which the initial As always had the same sign for 
all X. We have assumed however As = As, since the fit is 
unable to discriminate strange quarks from anti-quarks. 
This is really an assumption: unlike the spin-averaged 
case where the distributions s and s will be rather simi- 
lar (the integral of s — s has to vanish), there is a priori 
no need for As and As to have the same size or even the 
same sign. 

Qualitatively, the main features of our extracted 
strange sea distribution arise in the following way: the 
(kaon) SIDIS data, within the leading-twist framework 
we employ, turn out to prefer a small and likely positive 
As at medium x, while inclusive DIS and the constraints 
from /3-decays demand a negative integral of As and so 
force As to turn negative at low x. Given the importance 
of As, we address these constraints and their significance 
and implications in more detail in the following. 

We start by we analyzing the behavior of the trun- 
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FIG. 8: Same as Fig. |6] but for the truncated first moment 
of the polarized strange distribution As^'''' °°^~^ °l . 



cated first moment, As^^^" '^*'^"*^! , around the minimum 
defining the best fit. Figure [5] shows the increase of 
of the fit against variations of As^'t*^ *^"^^^] , along with 
the partial contributions of the various data sets. Evi- 
dently, the best fit has a truncated moment close to zero 
and only slightly negative, as we also saw in Tab. IIIII 
The shape of Ax^ around the minimum is dominated by 
the SIDIS data, and here primarily by the data for kaon 
production. All other data sets, pion SIDIS, inclusive 
DIS, and RHIC pp data, play less important roles, as 
expected (here one has of course to keep in mind that 
the impact of individual data sets seen in the Lagrange 
multiplier scans is always estimated in the "presence" 
of the other data sets, and therefore should not be con- 
strued as an independent fit result). As can be seen from 
Fig. [5] and Tab. IIIII the truncated moment of As remains 
close to zero if changes of Ax^ = 1 are permitted. For 
the more realistic choice Ax^/x^ = 2%, one finds that a 
much larger range of As^'^"""^^^' is allowed, extending 
from significantly negative to positive values. The size 
and even the sign of the considered truncated moment 
are, therefore, presently not well constrained. Nonethe- 
less, there is a trend for As(a;) to be positive at medium 
X ~ 0.1, even for the choice Ax^/x^ = 2% (see Fig. 2 
of ^). We note that the COMPASS experiment has re- 
cently presented a LO extraction of the polarized strange 
distribution from their kaon SIDIS data [79], which are 
not yet included in this work. These are consistent with 
small strangeness polarization down to below x ~ 0.01, 
but also allow a significantly negative As at a: ~ 0.005. 
Furthermore, an extraction of the integral of As over the 
range 0.02 < a; < 0.6 by the HERMES collaboration ^ 
yields 0.037 ± 0.019 ± 0.027, consistent with our result. 

We stress that beyond the "data-driven" uncertainties 
that we find for the polarized strangeness distribution, 
there could well be effects that are outside the leading- 
twist framework we are using here and that may have a 
significant impact on the extracted As. Given the gen- 
erally low hadron multiplicities in kaon events in the 
present SIDIS measurements, it is not ruled out that 
the kaon SIDIS data are affected by higher twist con- 
tributions and not suited for an extraction of leading 
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twist strangeness distributions. We note that the in- 
formation on the parton-to-kaon fragmentation functions 
in [s^l also primarily comes from unpolarized kaon SIDIS 
data and would not be reliable in the latter case either. 
A recent determination of the unpolarized strange dis- 
tribution in the nuclcon by HERMES from their SIDIS 
multiplicities shows an unexpected shape of the distribu- 
tion [80]. SIDIS measurements at smaller x, as well as at 
presently available x, but higher Q^, will likely be vital 
for clarifying these issues. These would become available 
at an electron-ion collider [sH - 

As can be seen from Fig. El the effects due to SU(2) and 
SU(3) flavor symmetry breaking in usage of the baryon 
semi-leptonic decay data, see Eqs. (|29)) . ([30)) . have only 
a very limited impact on the truncated moment of As. 
This, however, changes dramatically when the full first 
moment of As is considered, i.e., the contribution to its 
integral from x < 0.001. This region is presently not con- 
strained by any DIS or SIDIS data, but we remind the 
reader that the breaking parameters esu(2.3) come out 
very small, see Eq. (|34l) . in our analysis, as a result of 
the relatively small nominal uncertainty in the F, D val- 
ues, as we discussed in Sec. IIII Al This implies that the 
strange sea distribution will have a large and negative 
total first moment, AS^ — As^ -I- As^ — —0.114 as seen 
from Tab. IIVI which in turn can only occur if the distri- 
bution shows a sign change to negative values at small x, 
visible in Fig. [31 

It will clearly be important in the future to better un- 
derstand the strange contribution to nucleon spin struc- 
ture. If the full first moment AS^ is small, SU(3) symme- 
try in relating hyperon /3-decays to nucleon spin structure 
would have to be broken at the 40% level or so, which 
is not ruled out [63, H^. If, on the other hand, SU(3) 
symmetry is not broken significantly, the implication is 
that either As turns large and negative at small x, as in 
our fit, or that the present kaon SIDIS data do not allow 
a reliable extraction of As{x). On the theoretical side, 
there have been very recent lattice determinations of the 
integral AS^ [s^, which point to small values. Mod- 
els of nucleon structure, on the other hand, have led to 
quite varied predictions for the integral of As, ranging 
from small to large negative values |83|. We note that 
the "valence scenario" of [3l| has a first moment of the 
polarized strange distribution very close to zero, which is 
consequently at the expense of significant violation of the 
SU(3) flavor symmetry relation in Eq. (|30p. We finally 
stress that the size of As is not a topic of interest just for 
nucleon spin structure enthusiasts: as was pointed out 
recently [8J], the uncertainty in AS^ provides the single 
largest uncertainty in predictions of the spin-dependent 
elastic scattering cross sections of supersymmetric dark 
matter particles on protons and neutrons. 
Total quark and anti-quark spin contribution AS; In 
Fig. [51 we show the profile associated with variations of 
the truncated moment of the quark singlet distribution, 
^5^i,[o.ooi-.i] ^ J^^^^^dx[Au + Au + Ad+Ad + As + As], 
at = 10 GeV^. As expected, the main constraints 
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FIG. 9: Same as Fig. [6] but for the truncated first moment 
of the quark singlet distribution AE^'[° '"'^"^ °' . 



come from the DIS and SIDIS data. The value for the 
truncated first moment obtained in the best fit is signif- 
icantly higher than that for the full first moment given 
in Tab. lIV( which is a manifestation of the large negative 
contribution from strange quarks and anti-quarks that 
arises in our fit at small a;. Thus, keeping in mind the 
discussion about strangeness above, we conclude that if 
SU(3) flavor symmetry in relating hyperon /3-decays to 
nucleon spin structure is strongly broken, AS would be 
as large as ~ 0.36 or so, whereas it will be about 30% 
smaller if SU(3) holds well and the first strange moment 
ASs turns out to be large and negative. We note that 
such lower values of AS ~ 0.24 or so have usually been 
obtained in previous analyses relying on the use of SU(3) 
symmetry [31, 32, 33, 34, 36]. In any case, AS is certainly 
much smaller than the typical expectation of AS > 0.6 
in quark models. 

Spin- dependent gluon distribution Ag: We have already 
noted in our DSSV paper Q that the polarized gluon 
distribution Ag{x,Q^) comes out rather small in the 
presently accessed range of momentum fraction x, and 
prefers to have a node. At variance with the findings 
of Ref. [sil, we do not find any non-overlapping best- 
fit solutions with gluon polarizations of opposite signs. 
This duplicity is readily excluded by the RHIC pp data. 
The RHIC data in fact turn out to play a crucial role 
in constraining Ag [2^. The result is shown again in 
Fig. [3l We do not repeat the plot of the profile as a 
function of the truncated first moment of Ag here, which 
may be found in [28|. As can be seen from Tab. IIIIl 
the integral of A^ over the RHIC a:-region 0.05 to 0.2, 
^^RHic^ is found to be almost zero, while Tab. IIVI shows 
that extrapolation to all x results in the gluon spin con- 
tribution A5I = -0.084 at Q2 = 10 GeV^. We stress. 
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FIG. 10: Comparison of Ag/g for our best fit, at two repre- 
sentative , to the extracted Ap/g from photon-gluon fu- 
sion processes investigated by SMC [19], HERMES 18], and 
COMPASS HO, [in. These data were not included in our 
global analysis since a consistent NLO framework is not avail- 
able at present. 



however, that this result is not yet reUable due to the 
large uncertainty in extrapolation to x ^ 0. In any case, 
there are presently no indications of a sizable contribu- 
tion of gluon spins to the proton spin. This is in line 
with recent theoretical expectations obtained within an 
effective low-energy theory of broken scale invariance of 
QCD [855. Recent bag model estimates also point to rela- 
tively modest (but positive) values [s^ . Very large values 
of the integral of the spin-dependent gluon distribution, 
Ag^ ~ 1.5 or so at = 1 GeV^, as predicted based on 
considerations of the QCD axial anomaly [87| | , become in- 
creasingly disfavored, unless At; would show a steep rise 
at small x. Future data from RHIC for spin asymmetries 
in forward production of correlated hadron or jet pairs, 
and from running at 500 GeV c.m.s. energy, are expected 
to shed light on Ag at lower momentum fractions [ssj . 
Again, also a polarized electron-ion collider [8l| would 
be ideally suited to address this important question and 
to quantify the amount of gluon polarization at small x 
from measurements of scaling violations of the structure 
function gi. Other promising channels are, for instance, 
the polarizedphotoproduction of single- inclusive hadrons 
^ or jets [ia]. 

We have shown the comparison to some of the RHIC 
data in Fig. 0] (see also Ref. [11]). A way to access Ag 
in lepton-nucleon scattering is to measure final states 
that dominantly select the photon-gluon fusion process, 
heavy- flavor production, ip — > hX^ and tp h^h^ X, 
where the hadrons have large transverse momentum. 
Figure [TOj shows the corresponding results for the ex- 



tracted Ag/£ from SMC, HERMES, and COMPASS 
[isl [l^ . I20I [2l| , which have not yet been included in our 
global analysis. We also show in the figure our result, 
for two representative scales. It should be noted that 
this comparison is not quite consistent, as the extraction 
of Ag/g by the experiments was performed at LO level 
based on Monte-Carlo generators. Nonetheless, a small 
Ag at X ~ 0.08 — 0.2 as found in our analysis is also well 
consistent with the data from lepton-nucleon scattering. 
We expect that the data for the measured spin asymme- 
tries will be included in our global analysis in the future, 
after the NLO framework for them has been fully devel- 
oped and been compared to data for the corresponding 
spin-averaged cross sections. 



D. Exploring the fit parameter space 

In this Section we briefly present a few more details 
of the behavior of our total near its minimum, which 
has ramifications, in particular, for the use of the Hessian 
matrix method for estimating uncertainties. As we noted 
before, an advantage of the Hessian technique is that it 
allows to produce sets of "eigenvector PDFs" [11] , which 
in turn can be straightforwardly used in computations of 
other observables, in order to estimate their PDF uncer- 
tainty based on Eq. (jlip . For this, however, it is very 
important to know the range of validity of the method, 
i.e., to which degree is parabolic around its minimum. 
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FIG. 11: Correlations between the fit parameters {ui} and 
the eigenvector directions {zi}. The larger the box size the 
larger the overlap, see text. 



As described in Sec. IIIBl the first step in the Hessian 
method is to transform the fit parameters {ai} to a new 
set {zi} such that surfaces of constant turn into hyper- 
spheres in {zi} space, see Eqs. ((SJ, ^ [43]. Figure [TT] 
shows the overlap of each of the original fit parameters 



19 



{tti} with the eigenvector directions {zi}; the larger the 
box size the larger the contribution of a certain eigenvec- 
tor direction to a fit parameter a.;. The Zi are ordered 
in terms of the eigenvalues of the Hessian matrix: zi 
corresponds to the largest eigenvalue, i.e., a direction in 
parameter space where changes rapidly, whereas zig is 
only very weakly constrained by data. One can see that 
in many cases there is a fairly strong correlation between 
a given original fit parameter at and a single eigenvector 
direction z^. The parameters which appear to be con- 
strained best by current data are the normalizations of 
the sea quarks, Nu, N^, and Ng, esu(2,3) controlling the 
breaking of SU(2,3) symmetry, and au+u, C(d+d related 
to the small x behavior of Am -I- Am, Ad+Ad. Parameters 
determining the gluon distribution, Ng, ag, and rjg are 
less well constrained and mainly correlated with eigen- 
vector directions z-j to zi2. r/^, j^i^^, r/d+d^ ^^'^ Vs receive 
contributions from eigenvector directions which are only 
weakly constrained by data. As we shall see below, this 
general picture agrees rather well with results for the 
profiles for each fit parameter obtained with the Lagrange 
multiplier method. 

In Figure fT2l we investigate the behavior of around 
its minimum, making use of the transformed parameters 
{zi}. We vary one of the parameters Zi at a time, keeping 
all others fixed. Of course, since each Zi has in principle 
overlap with all fit parameters {a^}, the latter all vary 
in this procedure. The variation is done in such a way 
that a given change of Ax^ = T is produced. For truly 
quadratic behavior near the minimum, as is the under- 
lying assumption in the Hessian approach, the quantity 
— Axf, where Axf is the change in contributed 
by the parameter Zi that is varied, is trivially zero. This 
can be compared to the actual dependence of x^ on the 
varied parameter, making no use of the quadratic expan- 
sion in Any deviation of — Axf from zero will 
signal a departure from the quadratic behavior near the 
minimum. One can see from the figure that a choice 
Ax^ = 1 works reasonably well overall, in the sense that 
overall only fairly small deviations from zero occur. This 
implies that the Hessian matrix method is reliable for 
Ax^ = 1 and our eigenvector sets will produce faith- 
ful uncertainty estimates. Some eigenvector directions 
starting from zi2 and higher do show a certain depar- 
ture from the ideal behavior even for Ax^ < 1. This 
is most pronounced for zn to zig which are the least 
constrained parameters. In general we have found that 
the Hessian method breaks down rapidly once one goes 
beyond Ax^ — 1. Therefore we cannot provide eigenvec- 
tor sets corresponding to the more conservative error 
estimate Ax^/x"^ = 2% preferred in [28]. 

Figure [T5] shows the x^ profiles including the individ- 
ual contributions from the DIS, SIDIS, and RHIC pp data 
sets and from the F, D values for the fit parameters {a^}, 
obtained with the Lagrange multiplier approach. Clearly, 
while for some of the parameters the profiles are smooth 
and parabolic as expected in the simplest approach, for 
others they are not, showing not only non-parabolic be- 
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FIG. 12: Deviations from the expected parabolic behavior 
Ax^ = for the eigenvector directions {zi}, see text. Note 
that for better separation of the curves we have added an 
off-set i for each parameter Zi. 



havior but variously asymmetric shapes, multiple minima 
or almost flat regions. It is worth pointing out that these 
behaviors are not related to a lack of flexibility of the 
input parameterizations, but to features of the data it- 
self. For example, the double minima observed for iV^ 
and are associated with two possible "best-fit solu- 
tions" to the pion sidis asymmetries, which show strong 
fluctuations. 

In most cases, the behavior is still reasonably quadratic 
within Ax^ < 1, however, which further justifies the 
applicability of the Hessian method for Ax^ = 1. Be- 
yond that, simple extrapolation based on an assumed 
quadratic behavior may give misleading results. We re- 
call that the central values for the parameters can be 
found in Tab. [TTl 

Reducing the number of parameters of the fit would 
improve the constraints on the remaining ones, however 
at the expense of reducing the quality of the fit. The 
resulting constraint in that case would be strongly de- 
pendent on the functional form assumed for the PDFs. 
In this sense, the robust error analysis based on Lagrange 
multipliers allows to use more flexible functional forms. 
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FIG. 13: The profiles obtained with the Lagrange muhiplier approach for the parameters {ai} of the fit. The solid lines 
represent the increase in above the value obtained in the best fit. Dashed-dotted lines, long dashes, dots, and short dashed 
lines represent the partial contributions coming from inclusive, semi-inclusive, RHIC and baryon-decay data, respectively. The 
parameters /3i have been fixed here and are hence not included. 



IV. APPLICATIONS OF OUR UNCERTAINTY 
ESTIMATES 

The 38 eigenvector PDF sets that we have con- 
structed for the Hessian matrix method are ideally suited 
for estimating the PDF uncertainty of other observables. 
In this way, one can for example gauge the accuracy that 
future additional measurements will need to have in or- 
der to have a significant impact on our knowledge of the 
PDFs. In this Section, we will present a few examples 
for the case of RHIC. In view of the fact that RHIC has 



just completed its first physics run at V5 = 500 GeV, 
we will focus on predictions for this c.m.s. energy. 

Figure [H] shows the NLO double-spin asymmetries 
All defined in Eq. ((36|) for pp — > hX {h = n^,n^) and 
pp jetX, for our central DSSV fit (solid lines), includ- 
ing the Hessian uncertainty bands for Ax^ — 1 using 
Eq. pT|) . One can see that the asymmetry for 7r° re- 
mains very small until about px ~ 20 GeV, as could be 
expected from a simple scaling of the asymmetry All ai 
VS = 200 GeV shown in Fig. (gD with xt = 2pt/VS. 
It then rapidly increases. The asymmetry for negatively 
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charged pions remains small for all transverse momenta 
and in fact turns slightly negative at high p^- In con- 
trast, is higher than A^j^, reaching about 3% at 
the highest pr shown. The behavior of the various pion 
asymmetries is closely tied to that of the polarized gluon 
distribution: at high pT, relatively large values of x are 
relevant, where our Ay is positive, and quark-gluon scat- 
tering dominates. An important contribution to the spin- 
dependent cross section thus involves the combination 
{/^u®Dl + /^u®Dl + M®Dl + /^d®D'^)®l^g of par- 
ton distributions and fragmentation functions. For 7r+ 
production, the u quark and d anti-quark contributions 
are expected to dominate, as these are valence quarks in 
a 7r+. The combination Au -1- Ad is positive as Fig. [3] 
shows. The large negative contribution associated with 
Ad is suppressed here. For 7r° production, the participat- 
ing fragmentation functions are all equal, and one probes 
the sum of up and down quark and anti-quark distribu- 
tions, which is positive but smaller than Aw -t- Ad. Fi- 
nally, for vr^ production, the main contribution involves 
Ad-I- Aw, which explains the downturn of to negative 
values at high pT- Clearly, the three pion asymmetries 
are also sensitive to the sign of Ay. We note that pre- 
liminary results for at ^fS = 200 GeV have recently 
been reported from RHIC f27| . 

Similar features as for the tt'^ asymmetry are ob- 
served for jets. Very roughly, one finds that A^^\^{pt) ~ 
A'l^^{kpT), where k ~ 0.5 or so, corresponding to the fact 
that on average only the fraction k of the total jet mo- 
mentum is taken by an observed tt*^. This implies that, at 
a given px, the jet spin asymmetry is smaller than that 
for ttO. 

We have seen in the previous Section that the SIDIS 
data have given some first insights into the flavor struc- 
ture of the polarized sea distributions of the nucleon. On 
the other hand, the uncertainties in SIDIS are still quite 
large, and it is in particular difficult to quantify the sys- 
tematic uncertainty of the results related to the frag- 
mentation mechanism at the relatively modest energies 
available so far. Complementary and clean information 
on Alt, Am, Ad, and Ad will come from pp — > W^X 
at RHIC, where one will exploit the maximally parity- 
violating couplings of produced W bosons to left-handed 
quarks and right-handed anti-quarks [H, [9l| . The high 
scale set by the W boson mass makes it possible to ex- 
tract quark and anti-quark polarizations from inclusive 
lepton single-spin asymmetries in W boson production 
with minimal theoretical uncertainties, as higher order 
and sub-leading terms in t he p erturbative QCD expan- 
sion are suppressed 0, [ol, HJ [9^1 ■ 

As a further application of our Hessian uncertainty 
PDFs, we show in Fig. [13] the single-longitudinal spin 
asymmetries, 

Al^^4-^, (38) 

for the processes pp £^X, where the arrow denotes 
a longitudinally polarized proton and £ = e or ^. The 
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FIG. 14: Double-spin asymmetries All for jet and pion pro- 
duction at RHIC at ^/S — 500 GeV as functions of the trans- 
verse momentum pr of the produced final state. We show the 
results for the best-fit parton distributions from our global 
analysis, along with the uncertainties estimated using the Hes- 
sian method, allowing changes of one unit in x^- We also show 
the results for the "standard" scenario of ^l[ (dashed lines; in 
the lower plot the result for tt^ (tt") is given by the top (bot- 
tom) curve). We have used the CTEQ6M unpolarized parton 
distributions for the calculation of the denominator of the 
asymmetry. For the pions, we have assumed pseudo-rapidity 
coverage of \ri\ < 0.35, and for the jets of —1 < rj < 2. 



charged lepton is assumed to have been produced by a 
leptonic decay of the boson. The asymmetries are 
shown as functions of the charged lepton's rapidity 77iept, 
with 77iept counted positive in the forward direction of the 
polarized proton. We have integrated over px > 20 GeV, 
where px is the lepton transverse momentum. The re- 
sults shown in the figure are based on a simple LO cal- 



culation of the processes qq' 



the NLO 



corrections which we should in principle include for con- 
sistency are negligible for this observable [o^, [H, [13, [1^ ■ 
For production, neglecting all partonic processes 
but the dominant ud —^ one, the spin-dependent 
cross section in the numerator of the as ym metry is found 
to be proportional to the combination [93| 

A?2(a;i)d(a;2)(l-cos6')2-Ad(xi)?l(a;2)(l-l-cos6')2 , (39) 

where 6 is the polar angle of the electron in the par- 
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FIG. 15: Single-longitudinal spin asymmetries for charged- 
lepton production at RHIC through production and decay of 
W bosons. The bands correspond to our uncertainty esti- 
mates based on the Hessian Ax^ = 1 eigenvector PDFs. We 
also show in the figure the spin asymmetries obtained for the 
"standard" and "valence" scenarios of 31]. For the case of 
W~ , we also show the results of two fits for which the ratio 
Ad{x)/d{x) is forced to turn to -1-1 as a; ^ 1, see text. 



tonic c.m.s., with ^ = in the forward direction of 
the polarized parton. At large negative jyicpt, one has 
X2 S> xi and 6 3> 7r/2. In this case, the first term in 
Eq. (j39p strongly dominates, since the combination of 
parton distributions, Au{xi)d{x2), and the angular fac- 
tor, (1 — cos^)^, each dominate over their counterpart in 
the second term. Since the denominator oi is propor- 
tional to u{xi)d{x2){l — cos 6')^ -I- d{xi)u{x2)il + cosO)"^, 
the asymmetry provides a clean probe of Au{xi)/u{xi) 
at medium values of xi. Indeed, the Hessian uncertainty 
band for Am shown in Fig. ^ is directly reflected in the 
band we show in Fig. [151 We also show in the figure 
the spin asymmetries obtained for the "standard" and 
"valence" scenarios of [31]. The latter has a large and 
positive Au distribution at the relevant x ~ 0.1, which 
clearly shows in the asymmetry. By similar reasoning, at 
forward rapidity yyicpt ^ the second term in Eq. ([55]) 
dominates, giving access to —Ad{xi)/d{xi) at relatively 
high xi. We have discussed in Subsec. fill CI that there is 
interest in the question if the polarized down-quark dis- 
tribution turns positive in the large-x region, for which 



there are currently no indications, see Fig. ([5]). As Fig. [15] 
shows, the asymmetry for W~ production becomes large 
and positive at high ?7iopt, which precisely reflects the fact 
that Ad{x) remains negative at high x in our DSSV fit. 
It is interesting to investigate how the asymmetry might 
look if Ad/d were to turn to +1 as x ^ 1. In order to 
do this, we have produced two fits where Ad/d is forced 
to have this behavior. The two fits are characterized by 
the value xq where Ad{x,M^) changes sign from nega- 
tive to positive values. We have chosen xq = 0.67 and 
Xq = 0.55. The values for these two fits are of course 
significantly worse than for our DSSV best fit, by about 
four units for xq = 0.67 and about 25 units for xq = 0.55. 
The results for the two fits are shown by the dotted lines 
in Fig. [TSl It should be well possible at RHIC to mea- 
sure the asymmetry for values ryicpt out to > 2 [ssj . We 
note that the behavior of Ad/d at high x will also be 
further addressed by experiments at Jefferson Lab after 
the 12 GeV upgrade 96]. 

For W'^ production, one has the following structure of 
the spin-dependent cross section [93]: 

Ad{xi)u{x2){l + cos0f - Au{xi)d{x2)(l-cosef . (40) 

Here the distinction of the two contributions by consid- 
ering large negative or positive lepton rapidities is less 
clear-cut than in the case of W~ . For example, at nega- 
tive 77iopt the partonic combination d{xi)u{x2) will dom- 
inate, but at the same time 6* ^ 7r/2 so that the angular 
factor (1 -I- cosS)^ is small. Likewise, at positive ryiept 
the dominant partonic combination Au{xi)d{x2) is sup- 
pressed by the angular factor. So both terms in Eq. (liP)) 
will compete essentially for all ryicpt of interest. This is 
reflected in the behavior of the calculated spin asymme- 
try shown in Fig. [151 which does not show as clear 
features as the one for W~ bosons. Nonetheless, the W'^ 
measurements at RHIC will of course still be of great 
value. In fact, our global analysis technique is precisely 
suited for extracting information on the polarized PDFs 
even if there is no single dominant partonic subprocess. 

V. CONCLUSIONS 

We have presented details of a recent study of the 
helicity parton distribution functions of the nucleon, 
which used experimental information available from in- 
clusive and semi-inclusive polarized deep-inelastic lepton- 
nucleon scattering and from polarized proton-proton 
scattering at RHIC. The data sets were used jointly in a 
next-to-leading order global QCD analysis, which allows 
to extract the set of parton distributions that provides 
the optimal overall description of the data, along with 
estimates of its uncertainties. We have presented tech- 
niques and computational methods that speed up the 
next-to-leading order calculations for pp scattering to the 
level required in practice for a global analysis. Our tech- 
nique is formulated in Mellin moment space. A key fea- 
ture is that the computationally most challenging parts 
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are done only once, prior to the fit. Use of a Monte-Carlo 
sampling method allows us to perform this one-time cal- 
culation very efficiently. 

Our extracted parton distributions show particularly 
interesting features in the sea quark and gluon sector. 
We find evidence for a mostly positive Au and a negative 
Ac? distribution, so that Au — Ad is positive. This behav- 
ior has been predicted by a number of models of nucleon 
structure. The polarized strange quark distribution As 
comes out slightly positive at medium x, which is driven 
by the semi-inclusive kaon DIS data and could be subject 
to rather large systematic uncertainties. As turns neg- 
ative at cc < 0.02 as a result of constraints from SU(3) 
symmetry, which have a relatively small nominal error. 
If true, this means that As acquires its large negative in- 
tegral essentially completely from the small- a; region. As 
a further consequence, quark and anti-quark spins com- 
bined contribute about a fourth to a third of the proton's 
spin, with the lower value arising if strange quarks and 
anti-quarks are indeed strongly negatively polarized at 
low X. Finally, we have found that the gluon helicity dis- 
tribution Ag{x, Q^) is small in the region of momentum 
fraction accessed directly so far by RHIC, with likely a 
node and an almost vanishing integral over that region. 
Reliable statements about the full gluon spin contribu- 
tion to the proton spin are presently not yet possible. 

We have performed uncertainty estimates for our po- 
larized parton distributions, using both the Lagrange 
multiplier technique and the improved Hessian approach. 
To obtain these, a large number of additional fits are nec- 
essary, for which the computational techniques we have 
developed are particularly important. We find that both 
approaches yield consistent results for moderate depar- 
tures from the best fit, typically Ax^ — 1. For larger 
Ax^, significant differences develop as a result of depar- 
tures from parabolic behavior of around its minimum. 
This implies that the Hessian matrix method becomes 
unreliable. We have produced a set of 38 "eigenvec- 
tor" parton distributions for the Hessian method with 
Ax^ = 1 [93, which may be used to estimate the un- 
certainty of any observable that depends on the distri- 



butions. We stress, however, that we presently prefer a 
more conservative choice of Ax^/x^ = 2% as a tolerance 
criterion for acceptable parton distributions. Unfortu- 
nately, the behavior of around its minimum does not 
warrant use of the Hessian method for producing eigen- 
vector parton distributions in this case. 

We have used the Ax^ — 1 eigenvector distributions to 
obtain predictions for spin asymmetries for high trans- 
verse momentum pion and jet production in polarized 
proton-proton collisions at 500 GeV center-of-mass en- 
ergy at BNL-RHIC, as well as for W boson production. 
The former would give information on A^ at lower x, 
while the latter would provide a clean new probe of the 
polarized quark and anti-quark distributions, which is 
important in view of the uncertainties inherent in semi- 
inclusive DIS. Our results indicate that there is signifi- 
cant potential for RHIC to provide further important in- 
sights into nucleon helicity structure. It will be straight- 
forward to include all the forthcoming data in the global 
analysis. 
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